Dynamics of vortex glass phase in strongly type II superconductors 
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Dynamics of vortices in strongly type-II superconductors with strong disorder is investigated 
within the frustrated three-dimensional XY model. For two typical models in [Phys. Rev. Lett. 
91, 077002 (2003)] and [Phys. Rev. B 68, 220502(R) (2003)], a strong evidence for the finite 
temperature vortex glass transition in the unscreened limit is provided by performing large-scale 
dynamical simulations. The obtained correlation length exponents and the dynamic exponents in 
both models are different from each other and from those in the three-dimensional gauge glass 
model. In addition, a genuine continuous depinning transition is observed at zero temperature 
for both models. A scaling analysis for the thermal rounding of the depinning transition shows a 
non-Arrhenius type creep motion in the vortex glass phase, contrarily to the recent studies. 

PACS numbers: 74.25.Qt, 74.72.-h, 74.40.+k 
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I. INTRODUCTION 



The application of superconductors crucially depends 
on the high electric current density without dissipa- 
tion. However, the resistivity would always be nonzero 
even in the presence of pining centers. This conven- 
tional picture has been changed with the discovery of 
high-T c superconductors [l| and the progress in random 
field systems [2j. Similar to the spin-glass system, Fisher 
et al suggested that, for strong disorder, the system 
freezes into a genuine thermodynamic amorphous vortex 
glass CVG) phase with some kind of glassy long-range 
orders [3,0]. The VG phase in strongly type II supercon- 
ductors has attracted considerable attentions both exper- 
imentally and theoretically [f| during the past more than 
one decades. It is of practical significance that the VG 
phase is a true superconducting state with a vanishing 
linear resistivity by diverging energy barriers. On the 
fundamental side, it is closely related to an important 
class of phenomena in condensed matter physics, such as 
spin glasses, random field systems |2[, and charge-density 
waves in solids Q. 

The evidences to support the existence of a VG phase 
have been reported in many experiments by means of 
the dynamic scaling of the measured current-voltage (IV) 
dataQ- However, Strachan et al. have shown that a per- 
fect collapse of the IV data is not the sufficient evidence 
for a VG transition Q, since the critical temperature 
and the scaling exponents are not uniquely determined 
by this dynamic scaling. 

Theoretically, the XY gauge glass model||, HU EI 
has been extensively employed to study the VG tran- 
sition. The values for the critical exponents are sim- 
ilar to those obtained in some experiments [7|. How- 
ever, lacking some of properties and symmetries due to 
the absence of net magnetic fields, it is questioned to 
be a model of disordered superconductors in an applied 
filed [H, Q [H, QUE!]- Some realistic models have then 
been proposed recently, but the conclusions were quite 



contradictory. Continuous finite temperature VG transi- 
tions have been given by most models with various crit- 
ical exponents. It was also observed that the VG phase 
disappears if the screening of the vortex interaction due 
to the gauge field fluctuation is included [l8[. The sim- 
ulation of the London-Langevin model suggested no VG 
phase fl9j . 

Among all vortex models, the disordered three- 
dimensional (3D) XY model with net magnetic fields has 
provided both equilibrium and dynamical vortex phase 
diagrams in Type-II superconductors with weak disor- 
der [2(J, The low field (weak disorder) low tem- 
perature phase is in general regarded as a dislocation- 
free Bragg glass with a quasi-long-range order [20|, which 
was observed directly in a neutron experiment [22} . Sev- 
eral dynamical simulations on the vortex matter with 
rather low fields for weak disorder have been performed 
in this model [U H^, [24]]. By a anisotropic frustrated 
3D XY model with strong disorder in the coupling con- 
stants, Olsson [3] provided evidence for the VG tran- 
sition in the unscreened limit. The correlation length 
exponent v — 1.5 ± 0.3 was obtained, consistent with the 
3D gauge glass universality. Within an isotropic model 
with different choice of strong random-coupling distribu- 
tion, Kawamura fTol ] reported similar results for the VG 
transition independently. Although the obtained value 
v = 1.2±0.3is slightly smaller, which within the er- 
ror bar also suggests a common universality with the 3D 
gauge glass model. However, it was found later that a 
convincing scaling collapse for helicity modulus could not 
be achieved in Kawamura's model[16(, possibly due to the 
small effective randomness in the small system accessed. 
Furthermore, to the best of our knowledge, the dynam- 
ical study in the frustrated 3D XY model with strong 
disorder is so far lacking, which is however more relevant 
to experiments in the context of VG transitions. 

In this paper, based on resistively-shunted-junction dy- 
namics, we perform large scale dynamical simulations on 
the frustrated 3D XY models for two typical sets of pa- 



2 



rameters in Refs. [3, [HI- The glass transition tem- 
peratures and the critical exponents are estimated based 
on the simulated IV data. The depinning transition at 
zero-temperature and creep motion far below the glass 
transition temperature are also studied. The rest of the 
paper is organized as follows. Sec. II describes the models 
and dynamic method. In Sec. Ill and Sec. IV, the main 
results for the VG transition, the depining and creep mo- 
tion of vortices are presented, and some discussions are 
also carried out. Finally, a short summary is given in the 
last section. 



II. MODEL 

The frustrated 3D XY model on a simple cubic lattice 
is given bv [Hini 

H = Jij cos(0 4 - (f>.j - Aij), (1) 

where <pi specifies the phase of the superconducting or- 
der parameter on site i, — (27r/$o) f- A • dl with 
A the magnetic vector potential of a field B = V x A 
along the z axis, Jy represents the random coupling 
distribution. The average number of vortex lines per 
plaquette is denoted by / = l 2 B/&o, where / is the 
grid spacing in the xy plane and $ is the flux quan- 
tum. We choose two typical sets of parameters used 
by 01sson[3| and Kawamurafl5j|. For convenience, the 
models with these parameters are called Models I and 
II respectively. In Model I, the random pinning po- 
tential is introduced in the coupling strength in the xy 
plane Jy = Jo(l + pey), where ey's are independently 
Gaussian distributed with zero mean and unit variance, 
p represents the pinning strength. The coupling be- 
tween the xy planes is J z — Jo/r 2 , (T is the anisotropy 
constant). We typically choose p = 0.4 which models 
strong pinning strength, 1/r 2 = 1/40, and / = 1/5. 
Simulations of Model I are performed with system size 
L xy = 100, L z — 60 satisfying L xy /L z — 5/3, much too 
larger than those in Ref. (l4| . In Model II, the quenched 
randomness is put in the coupling constant Jy in all di- 
rections, which is distributed uniformly on the interval 
[0,2 J ]. The filling factor is chosen to be / = 1/4. The 
present simulations of Model II are performed with the 
system size L = 64 for all directions, considerably larger 
than those in Ref. (la ]. 

The Resistivity-Shunted- Junction dynamics is incorpo- 
rated in simulations, which can be described as 

crh ^-^ • • dH > 

— ^Jfa - 4>j) = -qt + j cxm - 2^ % > ( 2 ) 

3 1 j 

where J cx t,i is the external current which vanishes ex- 
cept for the boundary sites. The rjij is the thermal noise 
current with zero mean and a correlator (??y (i)?7y(i')) = 
2akBTS(t — t'). In the following, the units are taken of 
2e = Jq = h = a = = I- 



In the present simulation, a uniform external current 
I x along a;-direction is fed into the system, analogous to 
exerimentsQ- The fluctuating twist boundary condition 
[25l ] is applied in the xy plane to maintain the current, 
and the periodic boundary condition is employed in the 
z axis. In the xy plane, the supercurrent between sites i 
and j is now given by jffy — Jij sin (9i — 9j — Ay — ry • A) , 
with A = (A x ,A y ) the fluctuating twist variable and 
6i = <pi + Ti ■ A. The new phase angle 0i is periodic in 
both x- and y-directions. Dynamics of A a can be then 
written as 

^ a = Jji \- J i% + Vij] -I<x,a = x, y. (3) 

<ij>a 

The voltage drop is V — —LA X . 

The above equations can be solved efficiently by a 
pseudo-spectral algorithm [2l| due to the periodicity of 
phases in all directions. The time stepping is done us- 
ing a second-order Runge-Kutta scheme with At = 0.05. 
The equilibration of the simulation should be ensured be- 
fore the measurement. So most of our runs are typically 
(4 — 8) x 10 7 time steps and the latter half time steps are 
for the measurements. The detailed procedure in the sim- 
ulations was described in Ref. [2l| . Our results are based 
on one realization of disorder. The present system size is 
much too larger than those reported in literature, it is ex- 
pected to exist a good self-averaging effect. We have done 
two additional simulations with different realizations of 
disorder for further confirmations, and indeed observed 
quantitatively the same behavior. In addition, it is prac- 
tically difficult to perform any serious disorder averaging 
for such a rather large system. Actually, the results from 
dynamic simulations on 3D XY model in the recent litera- 
ture were also for a single disorder realization(2li I23I l24j . 
mainly due to the large system simulated. For the data 
points presented in the following figures, the statistical 
errors are smaller or comparable to the symbol sizes. 

III. VG PHASE TRANSITIONS 

First, we study the VG phase transition in these 
two models. In Model I, the VG transition tempera- 
ture T g is estimated to be 0.123 ± 0.008 in equilibrium 
simulations 14] . The IV characteristics is simulated at 
temperatures ranged from 0.08 to 0.15, which must cov- 
ers possible Tg. I n the equilibrium simulations of Model 
II, Kawamura[15j obtained T g = 0.81 by performing the 
finite-scaling analysis of the Binder ratio and the mean- 
square current. The similar simulations on Model II are 
performed at the temperatures ranged from 0.5 to 1.1, 
which also covers the possible T g . In simulations on both 
models, we try to probe the system at currents as low 
as possible for each temperature. Figs. 1(a) and 1(b) 
present the resistivity R — V/I as a function of current / 
at various temperatures for Models I and II, respectively. 
It is clear that, at lower temperatures, R tends to zero as 
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FIG. 1: Current-resistivity curves at various temperatures for 
(a) Model I and (b) Model II 



FIG. 2: Dynamic scaling of IV data at various temperatures 
for (a) Model I and (b) Model II 



the current decreases, which follows that there is a true 
superconducting phase with zero linear resistivity. While 
R tends to a finite value at higher temperatures, corre- 
sponding to an Ohmic resistivity in the vortex liquid. 
These observations provide a evidence of the existence of 
the VG phase in both models. 

Assuming that the vortex glass transition is continuous 
and characterized by the divergence of the characteristic 
length and time scales t ~ (z is the dynamic exponent) 
, Fisher, Fisher, and Huse [J] proposed the following dy- 
namic scaling ansatz to analyze the glass transition from 
a vortex liquid with ohmic resistance to a superconduct- 
ing glass state, 

TRS}- Z = *±(7£ 2 /T). (4) 

where £ oc| T/T g — 1 \~ u is the correlation length which 
diverges at the transition. ^ (x) is a scaling function with 
the + and - signs corresponding to T > T g and T < 
T g . Eq. (|1]) was often used to scale measured IV data 
experimentally [7| . 

Right at T g , the RI curve should show a power law 
behavior R oc I~ a where a = (z — l)/2, which provides 
a convexity-concavity criterion to identify the VG tran- 
sition temperature as well as the dynamic exponent z. 
As shown in Figs. 1(a) for Model I that the value of T g 
is in between (0.12 — 0.13), because in the low current 
regime the RI curves show convexity below T = 0.12 
and concavity above T — 0.13. The RI curves at other 
temperatures within (0.12,0.13) can be obtained by in- 
terpolations. The temperature at which the RI curve 
most close to the power law behavior is regarded as T g 
and the RI power law exponent at T g gives the dynamic 



exponent z. The error bars are estimated by obvious de- 
viation from the power law behavior. In this way, for 
Model I, we obtain T g = 0.124 ± 0.002 and z = 5.8 ± 0.3. 
The value of T g is consistent with that in equilibrium 
simulations [14j. By the similar method, for Model II, we 
get T g = 0.81 ± 0.01 , z = 2.5 ± 0.2. Interestingly, the 
present value of T g in Model II agrees well with that in 
equilibrium Monte Carlo simulations by Kawamura[15j. 

Once T g and z have been estimated, we can examine 
the IV data at different temperatures by the dynamical 
scaling. Fig. 2(a) shows that the data collapses well 
according to Eq. (0| if using the correlation length expo- 
nent v = 1.6 ± 0.1. The error bar is estimated by tuning 
the value of v until the collapse becomes poor evidently. 
The value of v is very close to v = 1.5 ± 0.3 obtained in 
Ref. [HI through equilibrium Monte Carlo simulations 
of Model I. Also as indicated in Fig. 2(b) that, using 
v = 1.2 ± 0.1, an excellent collapse according to Eq. (4) 
is achieved. The value of v also agrees quite well with 
v = 1.2 ± 0.3 estimated in an equilibrium Monte Carlo 
simulations of Model II [la]. Interestingly, although the 
values of v in both Models lie in the range [1.0 — 2.0] 
usually observed experimentally [3] , they are close to but 
different from each other. Since the present two models 
involve different symmetries (anisotropy) and different 
disorders included, in our opinion, it is not unlikely that 
they represent different universality classes. 

It should be mentioned that the present analysis 
method for the VG phase transition is not essentially 
inconsistent with that described in Ref. 8]. We also 
think that only the perfect collapse of the IV data is 
not a sufficient evidence for a VG transition, so we use 
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the convexity-concavity criterion to identify T g and de- 
termine z before performing the dynamic scaling. 

In equilibrium Monte Carlo simulations of Model II, 
some quantities failed to provide good scaling[l6|. The 
helicity modulus was used in the finite size scaling anal- 
ysis of the VG phase transitions in both models [Til fl^|. 
a nice scaling is obtained in Model I, but scaling fails ap- 
plied to Model II for data in system sizes L < 20. The col- 
lapse of the transverse helicity modulus with poor quality 
gives T g = 0.63^ = 1.5, which differed significantly from 
those in Ref. [la]. More seriously, it was impossible to 
collapse the data for the parallel helicity modulus. It has 
been observed [13] that the correct behavior required a 
great flexibility of the field induced vortex lines, which 
could be obtained either with a very large size or with 
weak interplane coupling along the field direction. For 
the isotropic system in Model II, the possible way to get 
a convincing scaling collapse of some quantities is to en- 
large the systems. In the present large scale dynamical 
simulations, a excellent collapse of the IV data in the 
dynamic scaling is indeed achieved. 

The exponents v in the present two models are close 
to v = 1.39 ± 0.20 evaluted by Olson and Young[ljl , 
but different from the recent more accurate result v = 
1.39 ±0.05 obtained by Katzgraber and Campbell [ll| in 
the 3D gauge glass model, suggesting that they are not in 
the same universality class. It follows that the difference 
in the quenched randomness and the introduction of net 
fields may change the static critical properties of the VG 
transitions. 

The dynamic exponents z in these two models are 
found to be quite different. The exponent z in Model 
I is high, in the range [4.0 — 6.0] usually measured 
in experiments 0, in Model II is however considerably 
low. Note that small values of the exponent z were also 
reported [la}- In addition, both exponents z in Models 
I and II can not fall even within the error bar of that 
in the 3D gauge glass model, which were estimated to be 
z = 4.2±0.6 in Ref. Q3 and z = 4.7±0.1 in Ref. [U, al- 
though the exponent z in Model I seems to be more close. 
It is possible that the disorder in the coupling constant 
along the filed direction in Model II reduces the effective 
pinning strength [l6| . leading to a small IV power law 
exponent at the VG transition. It is not expected that 
enlarging the system size further along the field direction 
would change the dynamic exponent z essentially. Nev- 
ertheless, the reason for the small value of z in Model II 
is not fully understood at the present stage, the further 
investigation is clearly called for. 



IV. DEPINING AND CREEP 

With the VG phase in hand, we then study the de- 
pining and creep motion of the vortices in this phase 
for both models. To study the depinning transition at 
zero temperature, we start from high currents with ran- 
dom initial phase configurations. The current is then 
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FIG. 3: Log-log plots of IV data at zero-temperatures for (a) 
Model I and (b) Model II. 



lowered step by step. The steady-state phase configura- 
tions obtained at higher currents are chosen to be the 
initial phase configurations of the lower currents in the 
next step. It becomes more difficult to measure the volt- 
age with the lower currents. In the vicinity of the critical 
current, a huge amount of the computer time is consumed 
to get accurate results. Fig. [3] exhibits the IV character- 
istics at T — for both models. Interestingly, we observe 
continuous depinning transitions with unique depinning 
currents [27]], which can be described as V oc (7 — I c )' 3 
with I c = 0.125 ± 0.001, j3 = 2.25 ± 0.02 for Model I 
and I c = 0.116 ± 0.002, /3 = 1.887 ± 0.01 for Model II. 
Note that the depinning exponents for both models arc 
greater than 1, consistent with the mean field studies on 
charge-density wave models [27j] . The depinning exponent 
is model dependent, possibly due to the different realiza- 
tions of the disorder. 

At low temperatures, the IV curves are rounded near 
the zero-temperature critical current due to thermal fluc- 
tuations. An obvious crossover between the depinning 
and creep motion can be observed around I c for both 
models at the lowest accessible temperatures. In order 
to address the thermal rounding of the depinning tran- 
sition, Fisher [271] first suggested to map this system to 
the ferromagnet in fields where the second-order phase 
transitions occur. This mapping was latter extended to 
the random-field Ising model[28] and flux lines in type-II 
superconductors (29j. If the voltage is identified as the or- 
der parameter, the current and temperature are identified 
as the inverse temperature and the field in the ferromag- 
netic system respectively, analogous to the second-order 
phase transitions, a scaling relation among the voltage, 
current, and temperature in the present model should 
satisfy the form 



V{T,I) = T 1 l 5 S[T- 1 / fl5 {l- I c /I% 



(5) 



where S{x) is a scaling function. 

It is implied that right at I = I c the voltage shows a 
power- law behavior V(T, I = I c ) oc T 1 / 5 and the criti- 
cal exponent 1/(5 can be determined. The log-log V — T 
curves are plotted in Fig. Q] (a) and (b) at three cur- 
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FIG. 4: Log-log plots of V — T at three currents around 7 c o 
for (a) Model I and (b) Model II. 



rents for Models I and II. In Fig. BJa), we can see that 
the critical current is between 0.115 and 0.135. The 
values of voltage at other currents within (0.115,0.135) 
can be evaluated by quadratic interpolations. The de- 
viation of voltage from the power law is calculated as 
the square deviations. The current at which the square 
deviation is minimum is defined as the critical current 
I c — 0.125 ± 0.02, consistent with those obtained at zero 
temperature. The temperature dependence of voltage at 
the critical current is also plotted in Fig. 0|a) . The slope 
of this curve yields 1/8 = 1.438 ±0.004. The similar anal- 
ysis in Fig. 4(b) yields I c = 0.116 ± 0.02 for Model II, 
consistent with that extracted from the zero-temperature 
simulation. The exponent 1/(5 = 1.227±0.003 is achieved 
by fitting the V — T curve in the low temperature regime 
at the critical current. 

With the critical exponent 8 and the critical current I c , 
we can adjust the depinning exponent (3 to achieve the 
best data collapse according to the scaling relation Eq. 
([5]) for I < I c . In Fig. [5] (a) and (b), a perfect collapse of 
the IV data at various temperatures below T g is shown 
with (3 = 2.25±0.01 for Model I and 1.89±0.01 for Model 
II. The values of (3 estimated from low temperature creep 
motion are in excellent agreement with those derived at 
T = depinning transition. Moreover, the scaling func- 
tion with the form V cx T 1 / 5 exp[A(l - ^)/T^ s } is de- 
rived in the creep regime for both models, which are also 
demonstrated in the legends of Figs. [5] (a) and (b). Note 
that the product of the two exponents (38 describes the 
temperature dependence of the creeping law. Interest- 
ingly, (36 ~ 1.56 for Model I and (36 ~ 1.54 for Model 
II are obtained, both deviate from unity, demonstrating 
that the creep law is a non-Arrhenius type. The values 
of (36 in both models are close to 3/2, which may moti- 
vate a further analytical work. In our opinion, it is not a 
coincidence that they are in the same universality class 
in the depinning transition. 

In a recent study of the depinning and creep motion 
of the flux line system in type-II superconductors [2{|, by 
simulations of overdamped London-Langevin model, Luo 
and Hu observed an Arrhenius law for the creep motion 
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FIG. 5: Scaling plot of the IV data at various temperatures 
below T g for (a) Model I and (b) Model II. 



with a linearly suppressed energy barrier for strong pin- 
ning |20j ] . inconsistent with the present study for strong 
disorder. It is worth noting that, in the London-Langevin 
model, the stable VG phase with the freezing of disor- 
dered vortex matter is not foundpjj and instead the vor- 
tices freeze like a window glass, called vortex molasses 
scenario. In the framework of the frustrated 3D XY 
model, the existence of a stable VG phase is well es- 
tablished in the unscreened limit in the present dynami- 
cal simulations, as well as in previous equilibrium Monte 
Carlo simulations [lil [HI [l6| . In the real strongly type-II 
superconductors, the screening induced rounding of the 
sharp VG transition is only a weak effect, and only visi- 
ble at temperatures very close to T g . The present good 
scaling behavior in the creep motion is just observed far 
below T g . We believe that the different nature of the 
phase in the London-Langevin model with strong pin- 
ning [29j and the VG phase in the present two models is 
the possible reason for the discrepancy. In addition, the 
Anderson-Kim creep law [32| is realized in the London- 
Langevin model with strong pinning {29l | , which may sug- 
gest that it is applicable to strong flux pinning in the 
conventional low T c superconductors rather than a VG 
phase with random point pins. 

The non-Arrhenius type creep behaviors have been also 
observed in charge-density waves with the mean-field re- 
sult (36 = 2/3 [30j| . the 3D random- field Ising model with 
136 « 3/2 Hi, (1+1) elastic interface with (38 « 2[Hj], and 
the flux line system in type-II superconductors for weak 
pinning in a Bragg glass phase with (38 ~ 3/2[2^|. It 
is surprising to note that the present combined exponent 
(38 w 3/2 in the frustrated 3D XY model for strong disor- 
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der is close to that in the 3D London-Langevin model for 
weak-pinning [2^ |. In the London-Langevin model of a 
fixed number of interacting particles, the vortex loop be- 
tween the planes perpendicular to the field is absolutely 
excluded, while in the present 3D XY models, the vortex 
loops between the planes can be induced by both ther- 
mal activations and the quenched disorder. So we argue 
that the disorder strengths in these two different kinds 
of models are hard to compare. Interestingly, the com- 
bined depinning exponent [36 ~ 3/2 was also observed 
in the depinning of domain walls in the 3D random-field 
Ising model [28|, possibly suggesting a universal rule in 
high dimensional elastic systems. Whereas the present 
results are different from recent results for (1+1) elas- 
tic interfaces in a disorder medium [3lj]. possibly owing to 
the two-dimensional nature in the latter. Further work 
is needed in order to clarify these observations. 

Note that the depinning of Bragg glass phase has been 
studied recently by Olsson using essentially the same 
model as Model I with a rather weak field / = 1/45 
[23| . The IV characteristics showed a unexpected behav- 
ior with a critical current that separates the creep region 
with an immeasurably low voltage at / < I c from a re- 
gion with V oc {I — Ic). This behavior is not observed 
in the present two models, possibly owing to the strong 
disorder and high fields. Their study together with the 
present one constitute a complementary picture for the 
depinning in the disordered 3D XY model with net fields. 

V. SUMMARY 

We have performed large scale dynamical simulations 
on the frustrated 3D XY models for strong disorder with 



two typical sets of parameters used in recent literature, 
within the resistively-shunted-junction dynamics. We 
first use the convexity-concavity criterion to identify T g 
and determine z, then perform the dynamic scaling on 
the simulated IV data. Adjusting the single parameter 
v, a perfect collapse of IV data is achieved for both cases, 
providing a evidence of the VG transition in the un- 
screened limit convincingly. Although the obtained cor- 
relation length exponents agree with the previous ones 
within error bars in equilibrium simulations, they are 
close to but different from each other, suggesting differ- 
ent universality class. New dynamic exponents are found 
to be parameter dependent. Both the static and dynamic 
exponents are different from the recent accurate results 
in the 3D gauge glass model. The nonlinear dynamical 
response far below the glass transition is studied system- 
atically. A non-Arrhenius type creep motion in the VG 
phase is observed, contrarily to the recent studies of the 
flux line system with strong pinning. The combined de- 
pinning exponent [38 = 3/2 is consistent with those in 
the 3D random-field Ising model and the flux line system 
with weak pinning, suggesting the common universality 
class in the depining transitions in these systems. 
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